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Strong field processes, which occur in intense laser fields are a test area for relativistic quantum 
field theory, but difficult to study experimentally and theoretically. Thus, modelling relativistic 
quantum dynamics and therewith the Dirac equation can help to understand quantum field the¬ 
ory. We develop a dynamic description of an effective Dirac theory in metamaterials, in which the 
wavefunction is modeled by the corresponding electric and magnetic field in the metamaterial. This 
electro-magnetic field can be probed in the experimental setup, which means that the wavefunc¬ 
tion of the effective theory is directly accessible by measurement. Our model is based on a plane 
wave expansion, which ravels the identification of Dirac spinors with single-frequency excitations of 
the electro-magnetic field in the metamaterial. We proof the validity of our relativistic quantum 
dynamics simulation by demonstrating the emergence of Zitterbewegung and verifying it with an 
analytic solution. 
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I. INTRODUCTION 


Fundamental properties in quantum electro-dynamics can be studied in strong field physics, in which for example 
intense laser beams can create electron positron pairs in the unstable vacuum[T]. The study of these processes 
analytically [2] or numerically [3HS] demands extended efforts. Therefore modelling systems are discussed in the 
literature, which allow to re-engineer physics in extreme conditions^ [7]. The capability of simulating intrinsic 
relativistic quantum processes like pair-creation or Zitterbewegung can therefore improve and test the understanding 
of quantum field theory. 

The latter mentioned Zitterbewegung, first considered by Schrodinger [8] , is a small quivering motion of particles in 
relativistic quantum mechanics. The theoretical existence of this quivering motion has been evidenced by numerical 
simulations of the Dirac equation and of quantum-field-theorv |5HTT| . However, the motion itself is very small and 
very fast. It can be estimated to be at the time scale of electron-positron pair creation seconds) and has 

an amplitude on the order of the Compton wavelength of an electron (10“^^ meters). Therefore, the experimental 
observation of this counter-intuitive property of nature is challenging. 

Nevertheless, the Zitterbewegung can be modeled in various other effective systems and it’s appearance has been 
subject of diverse investigations. It has been discussed for superconductors [T2j in 1970, for one-dimensional chains |T5] 
in 1990 and many considerations set in on semiconductors [THED] . In 2005 renewed discussions on Zitterbewegung 
in Sointronicspn 1^ . Graphene [2S1 |21j and carbon nanotubes PSI [55] raised up, (see also [57] 1. Zitterbewegung 
has also been investigated for trapped ions [58lj3T] . ultracold atoms|35] and Bose-Einstein condensates [33] with an 
experimental demonstration of the effect in 2010 [33]. More considerations on Zitterbewegung aroused for photonic 
crystals [35]. negative-zero-positive index materials [55] and binary waveguide arrays [37]. 

However, it is difficult to access the wavefunction itself in the listed experiments, where in particular higher order 
processes in quantum field theory are at least based on the knowledge of the wavefunction [35] . Recent experimental 
investigations demonstrated, that the wavefunction of the Dirac equation can be imitated by the electro-magnetic 
field of a waveguide structure with designable electro-dynamic properties [35]. Since this waveguide operates in the 
microwave regime, the simulated wavefunction can be directly probed in the experiment. But only quasi-stationary 
field configurations have been investigated and the question arises, whether the time-independent description in terms 
of frequency-eigensolutions is capable of forming dynamics, which correspond to the Dirac equation. 

Here, we extend the quasi-static theory to a dynamic description of an effective Dirac equation (section [IT]) and 
use it for the demonstration of the Zitterbewegung in theory. In the subsections jll Aj and jll Bj we repeat the already 
established description |39j of the Maxwell equations in metamaterials and show how solutions of the Dirac equation 
can be deduced. Next, we formulate a formal solution of the Maxwell equations by an expansion in plane-wave 
solutions in subsection m This solution can be identified with the unitary time-evolution of the Dirac equation, 
which is discussed in subsection jll D[ in which we also give an explicit mapping between the electro-magnetic field and 
the Dirac wave-function in frequency- and momentum space. The effective parameters for the mass and the speed of 
light of the emulated Dirac equation depend on the metamaterial an are derived in subsection jll E[ In section [Hi] we 
refer back to the Maxwell equations for deriving boundary conditions, with which electro-magnetic input pulses can 
be injected at the metamaterial interfaces in our simulation. 

After the introduction of the theory foundations, we describe the numerical implementation in section [rv] in which 
we also consider the properties of periodic boundary conditions, which are implied by our simulation method. 

In the Results section]^ we present the metamaterial simulations in three kinds of dynamical scenarios, in which 
we first consider the easiest possible setup for a Zitterbewegung with a Gaussian wavepacket excitation jV A| This 
setup is modified into a counterpropagating excitation, such that an experimental realization might be more feasible 


V B and finally we also account for the influence of the boundaries of the metamaterial V G 


The appendix contains a derivation of a formula of the expectation value of the position expectation operator, with 
which the Zitterbewegung can be computed semi-analyticalljj^ Another section discusses the absence of imaginary 
values in a real experiment, while the corresponding Dirac wavefunction consists of complex numberjB] 


II. THEORY DESCRIPTION 
A. The Maxwell equations 

The description of the one-dimensional Maxwell equations 

(la) 

(lb) 


-dxE^ = 

dxHy =-iuj€oer{uj)Ez . 
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in the photonic crystal is adapted from an effective model for the electro-magnetic field in a metamaterial waveguide, 
which has been developed recently [39] . Here, the electric and magnetic field is propagating in the x-direction and 
Eq and /io are the vacuum permittivity and vacuum permeability, respectively. The frequency dependent relative 
permittivity and relative permeability is given by 


er-(w) = — ( Co-^ ) and 

peo \ uj^Ld J 

^r{^) = ( Lq 2/^j ) ’ 

fj,o \ uj‘‘Cd) 


(2a) 

(2b) 


where we assume, that damping of the electro-magnetic field is negligible and the metamaterial is homogeneous along 
the x-direction in space. Damping terms may lead to a decay of electro-magnetic excitations, which might be subject 
to future studies. We assume specific material parameters for our model, which are listed in table |T] As a result, the 


TABLE I. Metamaterial parameters 

d = 8 mm element length 

p = 4 geometric factor 

C = 2.82 pF series capacitance of the loading elements 
Co = 58.8 pF/m per-unit-length capacitance of the trans¬ 
mission line segment 

L = 19.5 nH shunt inductance of the loading elements 
I/O = 314nH/m per-unit-length inductance of the trans¬ 
mission line segment 

This table contains a list of the parameter properties of the permittivity and permeability in equations §. 


electro-magnetic metamaterial properties assume the specific relations 

/ , /GHz\^ 

er{uj) = 1.66 — 188 | - | and 


= 1.00-141 




B. Effective Dirac equation 


(3a) 

(3b) 


As discussed in |39| , the Maxwell equations can be transformed to a one-dimensional, two-component Dirac equation 
by introducing 



Pi = 


(4a) 


P2 = 


(4b) 

If one identifies the effective mass 





'T 

II 

[Criu) - Prioj)] 

(5) 

and the effective energy 






[er(w) + Priuj)] . 

(6) 


one can rewrite the Maxwell equations, such that they are of the same structural form as the one-dimensional Dirac 
equation 


— icfxdx^p + m{uj)az^ = S{uj)(p . 


(7) 


ip = 


Here, ip is the two-component wavefunction 


P2 


( 8 ) 
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FIG. 1. (a) Dispersion relation of the metamaterial (111, in which the red and blue lines mark the upper band a;+ and lower 
band uj-, respectively. The grey-dashed lines are the functions f (w) and which are plotted along the a:-axis as functions 

of the frequency on the y-axis. (b) The same plot as in (a), but for the relativistic energy-momentum relation (A6l of the 
scaled, exact Dirac theory. The gray-dashed lines exhibit the scaled light cone |a;| = co\k\. 


c = l/^yeofJ-o is the vacuum speed light and and are the first and third Pauli matrices 


— 




(9) 


C. Time evolution of Maxwell equations 


In order to establish a dynamic Dirac theory, we first consider the time-evolution of the Maxwell equations ([^. 
Since the Maxwell equations are specified in frequency space, we evolve the electro-magnetic field in the metamaterial 
by the time-evolution of an expansion of plane wave momentum eigenfunctions Applying these plane waves 

yields the Maxwell equations in frequency- and momentum-space 


kE^^k{u}) = -w/roMr(w)i?j/,fc(aj), (10a) 

kHy^kiuj) = -ujeoer{u})E2^k{uj ), (10b) 

where the and Hy^ki'-^) are the electric and magnetic field amplitudes for a certain field frequency uj and a 

certain momentum k. From these Maxwell equations one can derive the dispersion relation 




( 11 ) 


The dispersion relation implies, that there exists a positive and negative momentum k for each frequency 
more, if one inserts the permittivity (2a) and permeability (pb|, one obtains a polynomial of second order 


means, that there exist two frequencies (upper band a;+ and lower band a;_) and the identical negative 
and —UJ- for each momentum k. This is plotted in figure (a) for the two positive frequencies uj+ and 
and lower band are divided by the band gap in which the product of er(a;)/rr(w) turns negative. Outside 
gap, the product €r{uj)^r{^) is always positive, and if also k^ is positive (ie. A: € K), it follows, that is 
and therefore w S K. The upper and lower boundary of the band gap in frequency space is given by the 
at which and fj,r are zero. This happens at the frequency 


ui. Further- 
in uj^. This 
values —a;_|_ 
u!-. Upper 
of the band 
positive too 
frequencies. 


wi = 10.4 GHz, 


at which is zero and at the frequency 


UJ 2 = 11.8 GHz, 


( 12 ) 


(13) 


at which is zero, for the parameters of table |T] The comparison of the dispersion relation of the metamaterial in 
Fig. 0 (a) with the relativistic energy-momentum relation of the scaled, exact Dirac theory in Fig. [^(b) illustrates a 
possible dynamical similarity of both descriptions around the momentum k = 0m“^. 
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Since the Maxwell equations are linear differential equations, the superposition of solutions is a solution again. 
Therefore, the time-evolution of the electro-magnetic field can be written down by the expansion 


E^{x,t) = 

E 

Ak, 


k 


Hy{X,t) = 

E 

Ak, 


k 

_} 



(14a) 


(14b) 


in which the sum runs over all discrete momenta k and all frequencies {uj+{k),uj-{k), —uj-{k)} for each 

momentum. The factor Ak is a measure for the summation in momentum space and is considered to be analogous to 
the measure dk of an integral. Since momentum space is spaced equidistantly in our approach, we write down a sum 
in (14) with measure Ak and focus on the plane wave solution in this more fundamental theory part. The spacing of 


Afc and it’s relation to position space is discussed in the section on numerical implementation IV 


equations (10), yielding 


In the expansion (14), the magnetic field can be expressed by the electric field by using one of the two Maxwell 




W^OMr(w) 


-Bz,fe(w), or 


OJeo€r{oj) 

— T ^z,k[^) ■ 


(15a) 

(15b) 


Therefore, the sum (14a) is already sufficient for a unique specification of the electro-magnetic field in the metamaterial. 


D. Time evolution of the effective Dirac equation 

In order to write down the time-evolution in Dirac theory, the plane wave solutions of the Dirac equation ([^ are 
to be considered. These are usually obtained by applying the plane waves at the Hamiltonian and solving the 
eigenvector and eigenvalue problem of the Hamiltonian matrix. The left-hand side of the Dirac equation ([^ turns 
into the matrix 


m{uj) k 
k —m{uj) 


(16) 


when applying the plane wave . The definition of the mass ([^ and energy ([^ together with the dispersion relation 
(11) ensure the identity 


£{uj)^ = k'^ + m{uiY , 


(17) 


which can be seen as relativistic energy-momentum relation of an electro-magnetic excitation in the metamaterial 
with wave number k and corresponding frequency w. Making use of 0 one obtains the eigenvector 


■(w) = 


^y2\£{uJ)\[\£{uJ)\ +m{uj)] 


with eigenvalue \£{io)\ and the eigenvector 

Wfe (w) = 


y^2\£{uj)\[\£{uj)\ + m{uj)] \\£iuj)\ + m{uj) 


\£{uj)\ + miuj) 
k 


-k 


(18a) 


(18b) 


with eigenvalue — (w)!, for the matrix (16). The solutions (18) are commonly known as bi-spinors of the Dirac equa¬ 
tion in one dimension [40) . We have normalized the bi-spinors in our notion, such that they form an an orthonormal 
base. Note, that the spinors (18) have the identical analytic expressions 


Wfe (‘^) = 


sign(fc) 


^y2\£{uj)\\\£{uj)\ - Tn{uj)] 


(19a) 
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and 


respectively, with the signum of k 


_ sign(fc) _ l'-\£{u})\ +m{uj) 

^ '' \/2|£(w)|[|f(a;)| - m{uj)] V ^ 


sign(fc) := 


1 ,iiO<k 
— 1 else. 


(19b) 


( 20 ) 


Even though the eigensolutions of the matrix (16) are determined by ( 18|19 ), one has to keep in mind, that this 
matrix is an integrated part of equation ([^. As a consequence the plane wave eigensolutions are formed by the 
electro-dynamics in the metamaterial. Since £{uj) is negative in the upper band and positive in the lower band, only 
the negative eigensolution appears in the upper band and only the positive eigensolution appears in 


the lower band. In fact, the spinors (18) can be expressed in terms of the electro-magnetic field. Keeping in mind, 
that £{oj) is positive at the lower band, the spinor can be written as 




1 ( ^i^-) + 

A/2£(w_)[£’(a;_) -I- m(a;_)] \ ^ 

_1_ ^-^r-(w-) 

y^[€r{uJ-) -I- IJ,r{u}-)]flr{uJ-) 


ck 

OJ — 


(21a) 


Similarly, £{u}) is negative at the upper band and by using the solution 11^ for Uj,, one can write 

sign(fc) _ f£{uj+)+ m{uj+)'' 

^2£’(a;+)[£:(w+) -h m{uj+)] V ^ 

sign(/c) 

ck I 5 
/ 


Wfc(w+)= 


^y[er{uJ+) + fj,r(,^+)]k-r{(^+) 


(21b) 


which is identical to (21a), except the sign of k. The sign(fc) factor in (21b) is necessary, because we want to keep 
the convention of the spinor solutions (18), as close to common relativistic quantum dynamics as possible (see also 
(A7)). Note, that the fraction of upper and lower component of the spinors ([M| corresponds to the factor (15a) 
between electric and magnetic field, if one accounts for the transformation law (13 of the electro-magnetic field of the 
Dirac wave function. Therefore, the free eigensolutions of the Dirac equation are linearly dependent to the plane wave 
solutions of the Maxwell equations and it just remains to determine the proportionality constant, to make a complete, 
unique link between the Maxwell equations of the electro-magnetic field in the metamaterial and the corresponding 
dynamic Dirac equation. To do so, we write down the time evolution of the wave function by expanding it with 


respect to the obtained eigenfunctions e®^“ and e 


— ^ikx 


of the Dirac equation ^ . 




+ ft (0.+) ft ft+) 


A/c, 


( 22 ) 


Here, (j)'^ (w) and <j)f, (w) are expansion coefficients, 


N = 


\(p{x, t)\'^dx 


(23) 


is a normalization constant and Afc is the measure for the momentum space summation, which is already mentioned 
in subsection IIC and further discussed in section IV Since (221 is a unitary transformation of the initial state of the 


wave function, the normalization constant will keep constant in time and <f{x, t) will always stay normalized at one. 
Note, that according to the considerations above, the positive eigensolutions which are proportional to will evolve 
with the lower band frequency a;_ and the negative eigensolutions which are proportional to u'jl will evolve with the 
upper band frequency uj+. 


If one compares the expansion of the electric field (14a) with the first component of the wave function (22) and 
requires equality for the prefactors of the exponential 6®®®“““+ one obtains 


£'z,fc(w+) = -(t>k (w+) 


sign{k)c^ 

•\/kr(w+) + /ir(a;+)| 


(24a) 
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Similarly, one obtains 


£;2,fe(w_) = (/)fc (w-) 


\/|e7(w^0'+”^r7(wI^ 


(24b) 


if one requires equality for the prefactors of 
with the second component of the wave function ( 
equivalent to equations (24), because the spinors 
the expansions (14). 


■. One may also equalize the expansion of the magnetic held (14b) 
2^ for the exponentials and . The result will be 

|2ip are linearly dependent on the electro-magnetic excitations in 


E. The metamaterial’s scaled Dirac equation 

According to the considerations in the sections |II C| and |II D| the metamaterial’s emulated quantum dynamics is 
expected to evolve as in Dirac theory, but natural constants like the electron mass and the effective speed of light 
depend on the metamaterial properties. 

For deducing these properties, we hrst consider the time-evolution equation 

ihdt(f{x, t) = Hip{x, t ), (25) 

with the Dirac Hamiltonian 

H = cpxCFx + moc^az , (26) 

in relativistic quantum mechanics. Here, mo is the electron rest mass and 

Px — ihdx (^^) 

is the momentum operator. After a Fourier transform of this equation into frequency space, the time-derivative idt is 
substituted by w, yielding 


Hip{x, cj) = hwp{x, cu). 


(28) 


If equation 0 could be cast in a similar shape with a linear uj at the right-hand side, it’s generalization to the time- 
domain would be straight forward and the mass and speed of light parameters of the emulated Dirac equation could 
be read of from the corresponding Hamiltonian on the left-hand side. This can be achieved by a Taylor expansion of 
the effective Energy £(uj) and effective mass m(uj) of the metamaterial Dirac equation ([^, at the frequency loq, at 
which £{uj) is zero. Thus, Wq is determined by setting the definition 0 oi £{uj) to zero and solving for w, resulting in 


Wo 


I 1 poC + p'^eoL 

CLd poCo + p'^eoLo ' 


(29) 


For the parameters chosen in table |T] wq has the value 11.0 GHz. The first relevant order of the Taylor expansion 
is sufficient for our consideration. For the mass m(w), the first relevant order is the zeroth order, whereas for the 
effective Energy £’(w), the first relevant order is the first order. The reason is, that Wq is chosen such that the zeroth 
order of £ (w) vanishes. Accordingly, one obtains 


- iaxdxP + m{uJo)(JzP = ^ 

for equation 0. If one defines the shifted frequency 

Aw := w — Wo , 


(w - Wo) (fi 


(30) 


(31) 


the effective speed of light 


m. \ 

9^ u, / 
^0 / 


CD ‘ = 


-1 


(32) 

















the new effective mass 


m' = m{u!o)- 


cd 


(33) 


and multiplies equation (301 with hco, then equation (301 can be transformed into 


PxCD'^x’f + in'c^azP = —hAu! . 


(34) 


This is consistent with the Hamiltonian in equation (28) and a negative time evolution, according to equation 
Comparing the corresponding Hamiltonian 


H = PxCDf^x + m'c^az 


(35) 


with the Hamiltonian of Dirac theory (26), one deduces the substitution 


C ^ Cjj ^ 
mo —?> m', 


(36a) 

(36b) 


which scales the Dirac theory of electrons such that it is similar to the emulated Dirac dynamics of the metamaterial. 

We point out, that equation (34) implies a negative time-evolution, because it’s right-hand side contains a minus 
sign in front of the hAunp expression. In a Dirac equation with positive time-evolution, this minus sign is a plus 
sign (see equation (28)). The minus sign originates from our chosen convention ( |4l^ . A different convention might 
require an inversion of the cc-axis. Therefore, if we want to compare the metamaterial dynamics with results from 
exact Dirac theory, we have to revert time. In particular we have to add a minus sign in the time argument of the 
semi-analytic expression (A16) and (A17) of the Zitterbewegung, which is derived in appendix [A| 


III. BOUNDARY CONDITIONS 

In this section, we consider electro-dynamic boundary conditions, for accounting for the finite space extension of 
the metamaterial. We assume, that the ’left’ end of the metamaterial is at location Xa and the ’right’ end of the 
metamaterial is at location Xh, where the notion ’left’ and ’right’ implies that Xa < Xb on the x-coordinate. As a 
consequence, the positions Xa and Xb are the limiters of three different regions, which are denoted by the indices (1), 
(2) and (3), in which the physical space is divided into 

X < Xa ■ region (1), the left input wave guide, 

Xa < X < Xb ■ region (2), the metamaterial, (37) 

Xb < X : region (3), the right input wave guide. 

The boundary conditions are obtained by integration of equation ([T ) over the infinitesimal interval [xa — e,Xa e] 
and [a;b — e^Xb + e]. We make the reasonable assumption, that the right-hand side of equation ([^ is bound and has a 
finite value. In this case the integrals over the intervals with infinitely small e imply, that the electric and magnetic 


fields must be smooth at the boundaries Xa and Xb, which means that 

E^^\xa,uj) = E^^\xa,u :), (38a) 

Ef'\xb,uj) = Ef\xb,L ^), (38b) 

H^^\xa,uj) = H^\xa,u :), (38c) 

{xb, uj) = {xb, uj) (38d) 


holds. 

We assume that the relative permittivity and permeability are just 1 for the propagation along the input wave 
guides of the regions (1) and (3), like it is for electro-magnetic waves in vacuum. Then, the dispersion relation 0 
simplifies to 


|w| = c\k\. 


(39) 
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For left propagating waves, in which = 

(141 can be simplified to 


-^2,left t) 


-^y,left t) 


II 

(i) 

0 and E^^_k{ 

-Uj) = Hy^_k{-Uj) 

= 0, the field expansion 

= E 

i?z.fc(a;)e*'=("- 

-ct) 

7 

(40a) 





= E 


-ct) 

5 

(40b) 

k,{LO,—Uj} 





for the regions (1) and (3). Similarly, for right propagating waves, in which Ez,-k{,(^) = Hy^-ki^) = 0 and Ez^ki—(^) = 
Hy k{—uj) = 0, the field expansion (14) can be simplified to 



(41a) 

k,{u),—uj} 


Hy,riM^,t)= Y Hy,k{uj)E'^^^+^*K 

(41b) 


k,{oJ, — Lo} 


for the regions (1) and (3). The simplified expansions imply a dispersionless translation of the signal 


^^ 2 ,ieft(ai, t) = E^^ieftix + cAt, t + At), (42a) 

^^z,right(ai, t) = E^^nght{x - cAt, t + At), (42b) 

-ffz,ieft(a:, t) = H^^ieit{x + cAt, t + At), (42c) 

-ff 2 ,right(a:, t) = i? 2 ,right(a; - cAt, t + At) (42d) 


in the regions (1) and (3). In other words: Once the left and right propagating input and output signals are known in 
the regions (1) and (3) at any position x, they can be trivially deduced from the above equations. It is most convenient 
to know the input signal at the metamaterial boundaries Xa and Xb. Therefore, we introduce the new coordinates 


The conditions (38) change into 


:= X — Xa for region (1) and 

(43a) 

:= X — Xb for region (3). 

(43b) 

Ei^\0,uj) = Ei^\xa,uj), 

(44a) 

Ei^\xb,u;) = Ei^H0,u;), 

(44b) 

3 

a 

II 

o' 

(44c) 

i/f)(x,,a;) = (0,cc) 

(44d) 


in terms of these new coordinates. If one inserts the field expansion (14) in these continuity conditions and keeps in 
mind, that the exponentials are linearly independent on the time interval ] — oo,oo[ for each frequency w, one 

obtains the boundary conditions in frequency space 


Ei%) + Ei%i^) = i?S(a.)e*'=- + 

(45a) 


(45b) 

<i(a.) + <1,(0;) = <i(o;)e*'=- + <l,(o;)e-'=- 

(45c) 

<i(o;)e^'=^‘ + <l,(o;)e-*'=“‘ = <i(a;) + <l,(o;). 

(45d) 


The relations (15) allow for the substitution of the magnetic field expansion coefficients in (45c) and (45d) with the 


electric field expansion coefficients, resulting in 

c^o 


CMo 


(46a) 
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physical boundaries 



input lines = 1, yUr = 1 

dump space 


FIG. 2. Implementation of boundary conditions within the periodic simulation area in subsection |V C| The unitary time- 
evolution according to (221 implicitly involves periodic bonndary conditions. Therefore, the physical boundary conditions 
according to section [ITT] are invoked at the space in-between the periodic boundaries of the whole simulation area. Since the 
simnlated metamaterial slab is embedded within some dump space, both boundary conditions are not interfering with each 
other. The left inpnt line is extended infinitely to the left and the right input line is extended infinitely to the right, respectively. 
They are described by the equations ( 40|41|42 1 of region (1) and (3). The green region of the physical metamaterial corresponds 
to region (2). 


Efl{> 










c/xo c/xq 


(3) 


where 


^(2) 


(w) 


k 

(2)7T 

W/Xo^f '{OJ) 


ujcoei^^ (lo) 
k 


(46b) 


(47) 


is an abbreviation for convenience. 


IV. NUMERICAL IMPLEMENTATION 


The numerical implementation makes use of an equidistant grid of states in momentum space, which by Fourier 
transform also implies an equidistant grid in position space. Assume, there are n grid points, each of them spaced 
by the extension Ax = 8 mm of the metamaterial’s unit element size (see table |T]). This implies a length I of the 
metamaterial, which is I = nAx. Since the plane waves should be 27r periodic after this extension I, the spacing 
in momentum space must be Ak = 2Tr/l. Thus, in momentum space, the momenta reach from —Ak{n — l)/2 
till Ak(n — l)/2, where n must be an odd number. The advantage of an equidistant momentum spacing is, that 


the plane wave spinors in the expansion (22) are complete and orthonormal basis functions on the spacial 


interval [—1/2,1/2], Furthermore, the simulation has periodic boundary conditions at the interval limits, because all 
exponentials are periodically continuing. Therefore, the simulations in the subsections | V A| and | V B| are valid, as long 
as the wavefunction’s probability density stays within the boundaries of the simulation area. 

In subsection I V C| the initial condition of the simulation is not specified at an initial time, but imposed by additional 
boundary conditions instead. These boundary conditions, which are discussed in section [Till are modeling the injection 
of input pulses at the physical boundaries of the metamaterials. However, the positions Xa and Xb, at which the 
injection boundary conditions are imposed are within the simulation area (see illustrative figure]^, such that they 
don’t interfere with each other. 


The time-evolution is computed according to equation (22), which is implemented numerically by first multiplying 
each expansion coefficient (j)k with it’s time-evolving phase factor at the certain time t and then summing over 

k. Since the sum over k runs over the exponentials and the discrete set of momenta k is equidistantly spaced, 
the sum is executed by performing a Fourier transformation of the (j>k array. 


V. RESULTS FROM SIMULATION 

In this results section, we apply the theory considerations from above, for demonstrating a Zitterbewegung of 
the emulated Dirac dynamics in the metamaterial. In a series of three different simulations, we first consider a 
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Gaussian wavepacket excitation as simplest possible setup in subsection V A For easier experimental implementation 
we also consider a moving Gaussian wavepacket in subsection V B In the last subsection V C we implement boundary 
conditions and account for the influence of these boundaries in the simulation. 


A. Gaussian wavepacket excitation 

1. Initial condition 

It is knownj^, that the Zitterbewegung only shows up for a simultaneous excitation of the positive eigen energy 
spectrum (above the mass gap of rriQC^) and the negative eigen energy spectrum (below the mass gap of moc^). 
Therefore, the simplest initial condition for a wavepacket, which exhibits Zitterbewegung dynamics is two Gaussian 
wavepackets in momentum space: One Gaussian wavepacket for the positive energy eigenstates and one Gaussian 
wavepacket for the negative energy eigenstates. Both Gaussians should be centered at momentum k = 0m“^. This 
description corresponds to the initial condition of the wavepacket 


4’k=^ I (48a) 

, (48b) 

at time t = 0 ns, where the positive- and negative energy eigenstates are excited equally. The width of the wavepacket 
should be of the order of the Compton wavelength h/{moc), which is the typical scale of the Zitterbewegung. By 
applying the metamaterial scaling (361, this turns into h/(m! cd)- Thus we choose the width in frequency space to be 


(Tfc = m'cuV^/fi. Correspondingly, the wavefunction (A5) takes the form 


ip{x) = 


VnJ-c 


77 / Akx I — Akx\ 

dk e + Uj^ e j e 


-(4)^ 


(49) 


with norm 


N = y/^ak ■ 


(50) 


2. Simulation 


The simulation is carried with a resolution of 401 grid points in momentum space. Therefore the simulation area 
has an extension of 3.2 m from the 1®* till the 401**' index, according to the numerical considerations section. 

We have plotted the probability density of the metamaterial simulation with initial condition (48) in figure in 
which we also compute the position expectation value 


/ <50 

dxx\ip{x,t)\'^ 

-OO 


(51) 


of the simulated probability density |(^(a:,t)p (red line). Equation (51) thereby shows the general formula for the 
position expectation, in which the limits of the infinite integration interval [— 00 , 00 ] are constrained down to the 
interval [—1.6 m, 1.6 m] in the case of the numerical simulation. For further investigation, we plot the red line of 
figure again in figure as solid black line. We compare this function with the position expectation value which 
is computed by using the simulated probability density of the exact Dirac equation (dotted line). The term ‘exact 


Dirac equation’ means, that the spinors and are replaced by the spinors and of equation (A7| and the 


frequencies a;+ and are replaced by -Isi£{k)/h of equation (A6) in the time evolution equation (22) of the Dirac 
equation. Another graph (plus-marked line) in figure]^ is the position expectation value, which is derived analytically 
in appendix]^ The final equation (A17) can be adapted to our problem by inserting the initial condition (48) and 
applying the scaling rule (|36|), which yields the semi-analytic equation 


{ip I x{t) ^ 




-2 W 


m 


sm — 


2f(^^ 


(52) 


Note, that we have inserted a minus sign in the argument of the sine function by hand, according to the negative 
time-evolution of the metamaterial simulation, which we have discussed at the end of subsection |II E[ 
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f/ns 


FIG. 3. Time evolution of a Gaussian wavepacket. The figure shows the probability density of the metamaterial simulation of 
the effective Dirac theory according to the dynamical evolution equation (221 with the initial condition (481, which is exhibiting 
a Zitterbewegung. Note, that the initial condition is specified at time t = 0 ns, which means that we show the computed 
backward time-evolution {t < 0 ns), being in one line with the forward time-evolution (0 ns < t). The red line displays the 
position expectation (511 and is also plotted as solid line in figure]^ 



simulation 
scaled Dirac 
analytic 


f/ns 


FIG. 4. Zitterbewegung of a metamaterial simulation. We show the position expectation value 


of the metamaterial 
The solid line is an 


simulation of the effective Dirac equation with initial condition ( |48[ ), which is illustrated in figure 
identical plot as the red line|^ Furthermore, the dotted line results from a simulation which is based on exact Dirac theory, as 
described in the main text. The plus marked line is a plot of the Zitterbewegung, given by the semi-analytic solution (52l. 


3. Discussion 

The line ‘scaled Dirac’ and the line ‘analytic’ in figure are both based on exact Dirac theory. For this reason, 
they are on top of each other and therewith identical. 

Furthermore, the comparison of the metamaterial simulation of the effective Dirac equation with the exact Dirac 
simulation in figureyields good agreement. The match of the ‘simulation’ line with the ‘scaled Dirac’ line and the 
‘analytic’ line is one of our main results. It implies, that our effective dynamic Dirac theory appears in metamaterial 
simulations such that the well-known Zitterbewegung can emerge. We conclude that metamaterials are capable of 
emulating the time-dependent Dirac equation in the case of well-suited parameters as for this setup. 


B. Counterpropagating excitation 

1. Initial condition 


The Gaussian wavepacket excitation in the above subsection is entering and exiting the simulation region very slowly, 
which means that the required electro-magnetic input pulse at the metamaterial interfaces will have an infinitely long 
head and tail. This is no useful property for an experimental realization and therefore, we demand a wavepacket which 
performs a Zitterbewegung but evolves more suitable in a prospective experiment. This property can be achieved by 
shifting the Gaussian wave packet (48) to the right in momentum space by the value fco = 20 m“^, implying the new 
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initial condition 




(pi =e \ -k j , (53a) 

_/ fc-fcp \ ^ 

= e V ‘"fc / , (53b) 

at time t = 0 ns. Here, the width of the wavepacket in frequency space is chosen to be Ufe = rn'cn/fi- Consequently, 
the wavefunction (491 changes into 




(54) 


where the norm (50) remains unchanged, except of course, the change of at- 


2. Simulation 


For the new setup, we increase the number of spacial grid points to 625. As a consequence, the simulation region 
is now extended over the interval [—2.496 m, 2.496 m]. 

Like in the above subsection |V A we have plotted the time-evolution of the probability density of the metamaterial 
simulation for the initial condition (53) in figure]^ One can see, that the positive and negative energy eigenstates are 
counterpropagating from x = —oo and a; = oo at time t = —oo, collide at a; = 0 at time t = 0 and move apart from 
each other to a; = — oo and a; = oo at time t = -|-oo. An oscillatory pattern is appearing at the colliding point (see also 
figure]^, which we attribute as Zitterbewegung dynamics. The reader might expect, that positive and negative states 
move in the same direction, because we have shifted the excitation by momentum fcp to the right in momentum space. 
But since the upper and lower band of the dispersion relation have the opposite slope at kg the shift in momentum 
space results in two counterpropagating wave packets. 

The metamaterial simulation of the effective Dirac equation and the simulation of the exact Dirac equation differ 
significantly from each other. Therefore, we plot them in two different plots (a) and (b) in figure]^ respectively. In the 
case of the effective Dirac theory of the metamaterial, there is an asymmetry which we explain with the asymmetry of 
the dispersion relation (see figure[^ (a)). The positive energy eigenstates are moving with group velocity {duj-/dk)\kg 
at ko of the lower band, while the negative energy eigenstates are moving with group velocity d{uj+/dk)\ko at ko of 
the upper band. Since the upper and lower band are differently shaped, the positive and negative energy-eigenstates 
are propagating at different speed through the medium. On the other hand, in the case of the exact Dirac equation 
in figure [^(b), the dynamics appears symmetric, which we explain with the symmetric upper and lower band of the 
relativistic energy-momentum relation in figure [^(b). As a result of the differences of the dispersion relations, there 
is an additional effective motion superimposed to the Zitterbewegung (see red line in figure]^ (a) and solid black line 
in figure]^ (a)), as compared to the wavepacket in the exact Dirac theory (see red line in figure]^ (b) and dotted black 
line in figure]^ (b)). 

Similarly as in the above subsection |V A[ we want to verify the exact Dirac theory by comparison with the semi- 
analytic solution in appendix [A| The new initial condition (53), substituted in equation (A17) yields 


/ I /.M \ 1 r -2(’^Y . f 2£{k)t\ 


(55) 


Note, that we inserted in a minus sign by hand in the argument of the sine function, as it has been done for equation 
(52) already. We plot the value of the solution (55) as plus-marked line in figure (b). 


3. Discussion 


The lines ^exact Dirac' of the exact Dirac simulation and the ^analytic' solution in figure [^(b) are coinciding. This 
is to expect, because the analytic solution is based on exact Dirac theory and is a validation of both descriptions. 

For a comparison of the metamaterial simulation with the exact Dirac theory, we subtract the value tj (49.7 ns) from 
the position expectation value in figure]^ (a), and plot it as the solid black line ‘modified simulation' in figure [^(b). 
We obtain very good agreement of both descriptions and take this match as further confirmation for the reliability of 
the simulation of the effective Dirac dynamics in the metamaterial. 
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FIG. 5. Time evolution of a colliding wavepacket. (a) The metamaterial time-evolution according to equation ( |22[ ) with initial 
condition (531 The red line marks the computed position expectation {x(t)) of ( |51[ ), which is also shown as solid line in hgure|^ 
(a), (b) The time-evolution of the same initial condition, but simulated with exact Dirac theory. The corresponding position 

expectation (red line) is also plotted as black dotted line in £gure|^(b). Like in figure]^ the initial condition for the simulation 
is specified at t = 0 ns, from which the backward and forward time-evolution is computed. 



FIG. 6. Zitterbewegung of colliding wavepackets. The solid black line in (a) and the dotted black line in (b), are the red lines 
in £gure[^(a) and respectively. The black cross-marked line shows the data from the semi-analytic formula of the position 
expectation value (55 I. For a comparison with the metamaterial simulated Dirac equation, we subtract the black solid line in 
(a) by t/ (49.7ns) and plot it together with the two other lines in (b). Once the steady motion is subtracted, the metamaterial 
simulation and the exact Dirac theory simulation agree with each other. 


C. Excitation input at the boundaries 

1. Initial condition 


In the case of a metamaterial with physical boundaries, the initial wavefunction (54) has to be replaced by an 


electric and magnetic field, which is propagating from the regions (1) and (3) into region (2) of the metamaterial. 
This can be done by specifying the incoming electric field of the regions (1) and (3) at the boundaries by 




(56a) 

(56b) 


with an arbitrary electric field amplitude, which is chosen to be if = 1 V/m, the frequency width cr^ = 0.52 GHz and 
the two frequencies uia = 13.81 GHz and uJb = 8.95 GHz. Here, the frequency uia is above the band gap and tOb is 
below the band gap, corresponding to the positive and negative excitation 


The expansion coefficients of the right propagating input pulse (56a) in region (1) are determined by the inverse 
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Fourier transform 




4 




(57) 


in time, with corresponding Fourier transform 

/ OO 

-OO 


(58) 


in frequency space. However, the equidistant grid of the simulated time evolution (22) in momentum space implies, 


that the sum of the Fourier transform has to be expressed in terms of an integral over the variable k instead of w. 
Such a momentum space integral of (58) in the fashion of (14) would read as 


fC 

^0 


dk 




(59) 


where the sum with momentum space measure Ak is replaced by an integral over dk. We consider a right propagating 
input pulse, therefore only right propagating modes with 0 < k are accounted for. Since we want to adapt the 
same density of states in region (1) as in region (2) in the numerical implementation, we relate frequency space to 
momentum space by the dispersion relation in the Fourier transform. Accordingly, the function a;(fc) implies an 
inner derivative in the integral 


rc 

E^P{x,t) = / 
Jo 


dk 


duj 

Ijk 




(60) 


and by comparison with (59), one obtains the relation 

du) 


E^Li^) = 


dk 


E'S(^) 




E'S^^) 


dk 

duj 




(61) 


By these relations, the expansion coefficients of the Fourier transform in time and frequency space, which is required for 
obtaining expansion coefficients of the input pulse (56a) can be related to the Fourier transform between momentum 


and position space (14), which is required for the metamaterial simulation and accordingly the Dirac time-evolution 

m. 


The same considerations also hold in an analogous way for the left propagating input pulse (56b) with expansion 
coefficients in region (3), as well as for the right propagating reflected pulse with expansion coefficients 

E^l(uj) in region (3) and the left propagating reflected pulse i?^^lj,(a;) in region (1). Note, that in the expressions 


(59) and (60) the integral over k is replaced by a sum over an equidistant momentum grid with momentum space 


spacing length Ak in the numerical implementation. 

Within the framework of the herein considered integral measures, the two reflected ‘output’ pulses and the electric 
field inside of the metamaterial are determined for each momentum, ie. each frequency according to the system 
of equations of the boundary conditions (45aI, (45b) and (46) of section III The fully determined electric, and 
therewith magnetic field in region (2) in form of the expansion coefficients E^^l and can be used to compute the 


wavefunction’s expansion coefficients (j)k by making use of the relations (24). Once the wavefunction is determined by 
this procedure, we normalize it by 


/ OO 

(w-)P + l^fc (w+)P] Ak, 

-OO 


(62) 


which according to Parseval’s theorem, is equivalent to normalization in position space (231. In contrast to the initial 
condition at time t = 0 along the whole x-axis in the above two subsections, this most realistic simulation specifies 


the ‘initial condition’ at the boundary positions Xa = —1.0m for (56a) and xt = 1.0 m for (56b I along the time-axis. 
However, through the relations ( |40|41|42| ), the alignment of the ‘initial 
the x-axis regions (1) and (3) at a certain time t <C 0 ns. 


condition’ can be recast from the time-axis to 
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FIG. 7. Time-evolution with boundary conditions. This figure shows a similar time-evolution as in 6gure[^(a), but this time, 
the wavefunction is given by the boundary electric heldj^at the two boundaries la = —1 m and = 1 m of the metamaterial. 
The wavepacket’s expansion coefficients and the time-evolution in region (2) are determined according to the procedure in 
subsection IID and section III The position expectation value (511 of this wavepacket is computed on the interval [xa,Xb] of 
region (2) and plotted as overlayed red line and in a magnified version in hgure[^ 


2. Simulation 


For the setup of this simulation we use 1001 grid points in position and momentum space, which implies that 
the simulation is extended on the interval [—4 m, 4m]. For this setup, the situation of figure applies, in which the 
boundary conditions of the metamaterial are implemented inside of the simulation interval, in which the physical 
metamaterial (corresponding to region (2)) is extended from Xa to Xb on the interval [—Im, Im]. The intervals 
[—4 m, — 1 m] and [1 m, 4 m] correspond to unphysical regions, which are depicted as dump space in figure ^ Therefore, 
even though the time-evolution has periodic boundary conditions, the wavepacket does not immediately reenter the 
simulation region (2) from the other side, after it went out before. 


3. Discussion 

The simulation of a Zitterbewegung with real metamaterial boundaries is shown in figure and the position 
expectation value of this simulation is overlayed as red line and plotted a second time in hgure|^ A similar simulation 
for the exact Dirac equation with the same initial condition differs much from the effective Dirac dynamics of the 
metamaterial (see also figure]^ and avoids the feasibility of a direct comparison. Therefore, there is no comparison 
with the corresponding dynamics of the exact Dirac equation in this subsection. 

In figure]^ one can see, that the electro-magnetic excitation in the setup is propagating from the physical boundaries 
towards the central region of the simulation area. In the overlap region of both excitations an oscillatory pattern 
emerges (see figure which we attribute as Zitterbewegung of the effective Dirac equation in the metamaterial. This 
demonstrates, that in the case of an ideal metamaterial, given by the equations Q and ([^ a Zitterbewegung of a 
simulated Dirac equation may occur for a real experiment. 


VI. CONCLUSION AND OUTLOOK 


The article discusses the feasibility on an exactly mappable, time-dependent Dirac simulation, which is actually 
emulated by electro-dynamics of Maxwell equations with designed electro-magnetic properties of a metamaterial 
structure. We present a time-evolution equation in a Dirac-like fashion (22) and give arguments of why it evolves very 
similar to exact Dirac theory. Furthermore, we consider scaling relations in subsection |II E[ with which the slower 
dynamics of the metamaterial can be mapped to the faster dynamics of elementary particles in Dirac theory. Based on 
the scaling rules and the exact Dirac equation, we have derived an explicit, semi-analytic expression for the position 
expectation of a general wavefunction in appendix [A| 

We demonstrate the dynamic description with three simulations: In subsection |V A| and |VB[ we compare our 
metamaterial simulations with exact Dirac theory and with the analytic solution, which is derived within this exact 
Dirac theory. We get very good agreement between the emulated Dirac dynamics and the exact Dirac dynamics, 
which we interpret as validation of our description. It is also a proof for the emergence of Zitterbewegung, which 
is a characteristic property of the free Dirac equation. Subsequently, we present a more realistic scenario, in which 
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FIG. 8. Zitterbewegung with boundary conditions. The plotted graph is the position expectation value, which is computed as 
red line in the boundary specified simulation in figure Here, we show a magnification of the region, in which the excitation 
of the positive energy-eigenstates coincides with the excitation of the negative energy-eigenstates in space. Like in figure an 
oscillatory pattern shows up, which we attribute as Zitterbewegung. 


boundary conditions are included, such that the simulation can be realized in experiment in subsection V C In this 
most realistic simulation we also observe an oscillatory pattern, which matches the Zitterbewegung. 

We conclude that a Dirac wavefunction can be emulated by the electro-magnetic field in a metamaterial, such that 
the Zitterbewegung occurs. Next steps are an experimental implementation of this theoretical study of the effect and 
also an investigation of the quantization of the system. 
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Appendix A: Zitterbewegung of a real electron in Dirac theory 


For verifying exact Dirac theory simulations, we consider the time-dependent position operator x(t) of the free Dirac 
equation in the Heisenberg picture in this section. This is already discussed in the literature. We simplify the already 
known derivation[Sl IHl I41| of the position operator for the one-dimensional case. Applying a general wavefunction to 
this position operator, yields a semi-analytic formula for the computation of the wavepacket’s position expectation 
value. 

We start out with the velocity operator ca^ in the Heisenberg picture 


2 

dtx = - [H, x] = ca^ , (Al) 

where H is the Hamiltionan of the exact Dirac theory ( |M| . A second application of the Heisenberg equation yields 
the acceleration operator 

dtcax = ^ [H, ctTa;] = - ^ - cTxH) . (A2) 

If one accounts for the time-independent constants H and p^, one finds the formal solution 

ccTxit) = {cax{Q) - c^PxH~^) . (A3) 

Integration of the velocity operator with respect to time gives the solution of the time-dependent position operator 
in the Heisenberg picture 


x{t) = Xq+ c^PxH - ^hc{axH ^-cpxH '^) (l - e ^ 

In the following, we deduce the expectation value of this operator for a general wavefunction 


(p{x,t) = 


1 




+ -+ ikx 


k— ^ikx 


) dk. 


(A4) 


OO 


(AS) 
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Note, that in contrast to (22) the spinors no longer have a frequency dependent mass in the Dirac theory of an 
elementary particle. Accordi 


ingly, due to the equations of motion (25) and (26) the energy-momentum relation 
£{k) = 


and the spinors 


ut = 


+ mQC^) ^ 


(fc)| -I- mo& 
' chk 


( —chk 


y2|£(fc)|(|f(fc)|+moc2) 


• TTIqC 


(A6) 

(A7a) 

(A7b) 


contain the natural constants mo, c and h. The computation of the expectation value of the position operator is as 
follows. First we define the function 


f{k) ;= (t)l u+ + (1)^ Ufc . 

With this abbreviation, the position expectation value of the position operator can be written as 


/ OO pOO pQ 

dx / dk f {k)^x{t) / 

-OO ^ —OO ^ —c 


dk'e’^'^fik'). 


(A8) 

(A9) 


The operator x{t) acts at on the right hand side integral, which is the integral over k'. This integral is a continuous 
sum over the exponential functions “ and one can commute the operator with the right hand side integration, 
where it solely acts at these exponentials. Hence, the position operator (A4) turns into a function of k', which means, 
that all momentum operators Px of the operator x{t) and H in x{t) turn into the number hk'. It remains a triple 
integral over the function 

k')f{k') (AlO) 

with the three integration variables a;, k and k'. If one performs the integration over x first, the exponential 

turns into the delta function 2tt 8{k — k'). Another integration over k' gives contributions only at fc = k', such that 

equation (A9) turns into 



poo 

{if\x{t)\ip) = 271 / dk f{k)^x{t,k)f(k'). 

J —OO 

(All) 

The remaining integrand can be reshaped by making use of the eigenvalue relations 



H{k)ul= £{k)ul, 

(AI2) 


H{k)u- = -S{k)ul , 

(AI3) 

the orthonormal property 

: 

(Al4a) 



(AI4b) 



(AI4c) 


%= 1 

(AI4d) 

and the basic identities 

~ 4 - chk 

Ut, <7iuZ = . 

^ ^ £{k) 

(AI5a) 


-+t 

% -1% = , 

(AI5b) 


--t ~+ 

= m ’ 

t chk 

(AI5c) 

(AI5d) 

Applying these identities, canceling 

some terms and rearranging the remaining terms yields 
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1 ( r^hhf 

{ip I x{t) \v’) = j dk}^Xo + [|(/'+(fc)P - \(l)~{k)\^] 


^ i mohc^ ' 


<^+(fc)*(/>-(fc) (l - - (!)+{k)(j)-{k)* (l - |. (A16) 


2 £{kY 

With a further usage of Euler’s formula this result can be further rewritten in the form 


{ip\x{t)\ip) = ^ j dkiyXo+[\(j)+{k)\^-\(j) (fc)p] 

mohc^ 


^hkt 


£{kf 


Im [(j)'^{k)*(j) (fc)] 


m 

(2£{k)t\ 


— Re (k)*(j) (fc)] sin ' ’ 


\ ^ 


)J} 


(A17) 


Appendix B: The wavefunction’s imaginary part in reality 

The whole description of the simulation of the electro-magnetic waves in the metamaterial and even the basic 
equations Q , on which our theory is built on, as well as in the description of the Dirac equation makes use of complex 
numbers. As a result, the electric and magnetic field has an imaginary part. However, in reality nothing is imaginary 
and there must be an answer of what happens with this kind of ‘degree of freedom’, which is not explicitly showing 
up in the real experiment. 

Consider the expansion 


(Bl) 

k 

The expansion will be exclusively real, if all expansion coefficients only have a symmetric real part 

Re (f{k)) = Re (f{-k)) (B2) 

and an anti-symmetric imaginary part 

Im (^f{k)) = -Im (^fi-k)) ■ (B3) 

On the other hand, the expansion will be exclusively imaginary, if all expansion coefficients only have an anti-symmetric 
real part 


Re (/(fc)) = -Re (/(-fc)) 


(B4) 


and a symmetric imaginary part 


Im (^f{k)^ = Im (/(-fc)) . 


(B5) 


In (14) the sum in the expansion runs over terms with negative frequencies and negative momenta at the same time 
and the discussion of exclusively real or imaginary functions in analogy to (Bl) gets more involved. Table |b] lists 


the different combination of symmetries, which are possible for the expansion coefficients ^^nd Hy k- Note, that 
the real part and the imaginary part can be chosen independently from each other. This means, that the resulting 
function will be only exclusively real or exclusively imaginary, if the symmetries of the real part and the symmetries 
of the imaginary part are both chosen either real or imaginary at the same time. According to these considerations, 
it is possible to choose expansion coefficients, such that the approximated function is exclusivel y re al. 

Having said that the electric and magnetic field can be exclusively real in the expansion (Bl), we additionally 


point out, that the Maxwell equations in momentum- and frequency space (10) are flipping the symmetry type from 


symmetric to anti-symmetric and vice versa with respect to frequency and momentum simultaneously, when they 
are coupling the electric expansion coefficients E^^k to the magnetic expansion coefficients Hy k- By comparing this 
flipping operation with table one realizes, that functions which are real, are only mapped to functions which are 
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TABLE II. Expansion coefficient symmetries 


coefficient 

frequency 

momentum 

function 

real part 

symmetric 

symmetric 

real 

real part 

symmetric 

anti-symmetric 

imaginary 

real part 

anti-symmetric 

symmetric 

imaginary 

real part 

anti-symmetric 

anti-symmetric 

real 

imaginary part 

symmetric 

symmetric 

imaginary 

imaginary part 

symmetric 

anti-symmetric 

real 

imaginary part 

anti-symmetric 

symmetric 

real 

imaginary part 

anti-symmetric 

anti-symmetric 

imaginary 


This table lists the different symmetries of the real part and the imaginary part of the expansion coefficients and Hy^k in 
the sum ( |14[) w ith respect to frequency uj and momentum k, in analogous extension to the symmetry considerations of the 
expansion ( |B1[ ). The first column ‘coefficient tells if the real part or the imaginary part of the wave function is considered, 
the second ‘frequency’ and third ‘momentum’ column tells, whether the real/imaginary part of the expansion coefficients are 
assumed to be fully symmetric or fully anti-symmetric with respect to frequency or momentum, respectively. The last line 
‘function’ tells if the resulting function of this expansion will be exclusively real or exclusively imaginary. 


real again. As a consequence, it is possible to choose expansion coefficients, such that the electric and magnetic helds 
are real in space and time and furthermore are compatible with the Maxwell equations (101. So the theory, which we 


have developed in the subsections above is compatible with a real electric held and a real magnetic held, even though 
complex numbers are showing up in the equations. 

However, the Dirac equation ( 25|26 ) explicitly implies complex wavefunctions. The time evolution (22) of the 
effective Dirac equation can be directly translated to the electric and magnetic held by the relations ([^ in real space 
and (24) in momentum and frequency space. Therefore, the question why the electro-magnetic held does not contain 


an imaginary part in reality, but in contrast the Dirac wavefunction does contain an imaginary part is still to be 
answered! 

The difference of the description of the Maxwell equations and the Dirac equation can be found in the expansions 


(14) and ( [22| ). While the sum in (22) only runs over the positive frequencies w_|_ and a;_, it also runs over the negative 
frequencies —a;+ and —ui- in (14). Therefore, on one hand, the simulated Dirac equation only makes use of the 


upper two bands of the dispersion relation © and is consistent with the two bands of the positive and negative 
energy eigenstates of the effective Dirac equation ( [M| ) (see also hgure . On the other hand, the dispersion relation 
(11) has four bands, of which the two negative ones have the just the negative value of the positive ones. So the 


expansion of the electro-magnetic field runs additionally over these two negative energy bands. As a conclusion, the 
time-evolution of the Dirac equation (22) omits the sum over the negative energy bands of the electro-magnetic field, 
without explicitly mentioning it. This omittance is introduced to be consistent with the two bands of conventional 
Dirac theory. However, the real- and imaginary part of the expansion coefficients fe(—w), Hy k{—oj) of the negative 
energy band of the Maxwell equations can be chosen according to the symmetry scheme of table [TJ such that the 
electric and magnetic held, which is transformed from the wavefunction by the relations (24) becomes exclusively 


real. If one eliminates the imaginary part by the above discussed symmetry procedure and divides by a factor of two, 
one acted as if one just took the real part of the wavefunction. Therefore, the considerations of this appendix are 
justifying that the real part of the simulated Dirac wavefunction is proportional to the electric and magnetic held of 
the underlying Maxwell equations, which can be measured in experiment. 
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